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Abstract 


This  paper  presents  a  multi-layered/multi-director  and  shear-deformable  finite  ele¬ 
ment  formulation  of  shells  for  the  analysis  of  composite  laminates.  The  displacement  field 
is  assumed  continuous  across  the  finite  element  layers  through  the  composite  thickness.  The 
rotation  field  is,  however,  layer-wise  continuous  and  is  assumed  discontinuous  across  these 
layers.  This  kinematic  hypothesis  results  in  independent  shear  deformation  of  the  director 
associated  with  each  individual  layer  and  thus  allows  the  warping  of  the  composite  cross- 
section.  The  resulting  strain  field  is  discontinuous  across  the  different  material  sets,  thereby 
creating  the  provision  that  the  inter-laminar  transverse  stresses  computed  from  the  constitu¬ 
tive  equations  can  be  continuous.  Numerical  results  are  presented  to  show  the  performance 
of  the  method. 
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1.  Introduction 

In  the  last  two  decades,  composites  have  found  increasing  application  in  many  engineer¬ 
ing  structures.  Recent  advances  in  the  technologies  of  manufacturing  and  materials  have 
enhanced  the  current  application  of  composite  materials  from  being  used  as  secondary  struc¬ 
tural  elements  to  becoming  primary  load-carrying  structural  components.  Due  to  the  inherent 
inhomogeneity  and  anisotropy  of  the  materials,  analysis  of  these  composite  structures  imposes 
new  challenges  on  engineers  [10,  11]. 

Plate  and  shell  structures  made  of  laminated  composite  materials  have  often  been  modeled 
as  an  equivalent  single  layer  using  classical  laminate  theory  (C.L.T)  [5,  8,  12,  13,  16]  in  which 
the  thickness  stress  components  are  ignored.  C.L.T.  is  a  direct  extension  of  classical  plate 
theory  in  which  the  well  known  Kirchhoff-Love  kinematic  hypothesis  is  assumed  enforced. 
This  theory  is  adequate  when  the  thickness  (to  side  or  radius  ratio)  is  small.  However,  lami¬ 
nated  plates  and  shells  made  of  advanced  filamentary  composite  materials  are  susceptible  to 
thickness  effects,  because  their  effective  transverse  moduli  are  significantly  small  compared 
to  the  effective  elastic  modulus  along  the  fiber  direction.  Furthermore,  the  classical  theory  of 
plates,  which  assumes  that  the  normals  to  the  mid-plane  before  deformation  remain  straight 
and  normal  to  the  plane  after  deformation,  under  predicts  deflections  and  over  predicts  natu¬ 
ral  frequencies  and  buckling  loads.  These  discrepancies  arise  due  to  the  neglect  of  transverse 
shear  strains.  The  range  of  applicability  of  the  C.L.T.  solution  has  been  well  established  for 
laminated  fiat  plates  [10,  12,  20].  In  order  to  overcome  the  deficiencies  in  C.L.T.,  refined 
laminate  theories  have  been  proposed  [3,  10,  11,  18,  19,  20,  21,  24].  These  are  single  layer 
theories  in  which  the  transverse  shear  stresses  are  taken  into  account.  They  provide  improved 
global  response  estimates  for  defiections,  vibration  frequencies  and  buckling  loads  of  mod¬ 
erately  thick  composites  when  compared  to  the  classical  laminate  theory.  A  Mindlin  type 
first-order  transverse  shear  deformation  theory  (S.D.T.)  was  first  developed  by  Whitney  and 
Pagano  [24]  for  multi-layered  anisotropic  plates,  and  by  Dong  and  Tso  [6]  for  multi-layered 
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anisotropic  shells.  Both  approaches  (C.L.T.  and  S.D.T.)  consider  all  layers  as  one  equivalent 

single  anisotropic  layer,  thiis  they  can  not  model  the  warping  of  cross-sections,  that  is,  the 
in-plane  distortion  of  the  deformed  normal  due  to  transverse  shear  stresses.  Furthermore, 
the  assumption  of  a  non-deformable  normal  results  in  incompatible  shearing  stresses  between 
adjacent  layers.  The  later  approach  also  requires  the  introduction  of  an  arbitrary  shear  cor¬ 
rection  factor  which  depends  on  the  lamination  parameters  for  obtaining  accurate  results. 
It  is  well  established  that  such  a  theory  is  adequate  to  predict  only  the  gross  behavior  of 
laminates. 

The  exact  analyses  performed  by  Pagano  [12,  13,  16]  on  the  composite  fiat  plates  have 
indicated  that  the  in-plane  distortion  of  the  deformed  normal  depends  not  only  on  the  lam¬ 
inate  thickness,  but  also  on  the  orientation  and  the  degree  of  orthotropy  of  the  individual 
layers.  Therefore  the  hypothesis  of  non-deformable  normals,  while  acceptable  for  isotropic 
plates  and  shells,  is  often  quite  unacceptable  for  multi-layered  anisotropic  plates  and  shells 
that  have  a  large  ratio  of  Young’s  modulus  to  shear  modulus,  even  if  they  are  relatively 
thin  [4,  18,  19,  20].  Thus  a  transverse  shear  deformation  theory  which  also  accounts  for  the 
warping  of  the  deformed  normal  is  required  for  accurate  prediction  of  the  elastic  behavior 
(defiections,  thickness  distribution  of  the  in-plane  displacements,  natural  frequencies,  etc.)  of 
multi-layered  anisotropic  plates  and  shells. 

In  view  of  these  issues  a  variationally  sound  theory  that  accounts  for  the  3-D  effects,  allows 
thickness  variation,  and  permits  the  warping  of  the  deformed  normal,  is  required  for  a  refined 
analysis  of  thick  and  thin  composites.  In  this  paper  we  have  endeavored  to  address  these 
issues  and  develop  a  new  finite  element  formulation  for  laminated  composites.  The  assumed 
displacement  field  is  continuous  through  the  thickness,  while  the  rotation  field  is  layer-wise 
continuous  (in  2-D)  and  can  be  discontinuous  across  the  finite  element  layers.  This  dis¬ 
placement  field  fulfills  o  priori  the  geometric  continuity  conditions  between  contiguous  layers. 
Furthermore,  the  assumed  displacement  field  is  capable  of  modeling  the  in-plane  distortion  of 
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the  deformed  normal,  without  increasing  the  order  of  the  partial  differential  equations  with 

respect  to  the  first-order  transverse  shear  deformation  theory.  In  this  formulation,  at  most, 
only  first  derivatives  of  displacement  and  rotation  fields  appear  in  the  variational  equations. 
The  practical  consequence  of  this  fact  is  that  only  continuity  of  finite  element  functions  is 
required  which  is  readily  satisfied  by  the  Lagrange  family  of  elements  [7].  Because  of  the  3-D 
features  in  the  formulation,  it  can  accurately  model  the  inter-laminar  conditions  and  predict 
the  3-D  edge  effects.  Finally,  like  the  single-layer  shear  deformable  theories,  the  proposed 
formulation  provides  fiexibility  in  the  specification  of  the  boundary  conditions  [7]. 

A  summary  of  the  paper  is  as  follows.  Section  2  presents  the  assiunptions  inherent  in  the 
proposed  model.  The  geometric  and  kinematic  hypothesis  of  the  multi-layered  shells  that 
accommodate,  to  an  appropriate  degree  of  approximation,  the  effects  of  transversal  warping 
of  the  cross-section  due  to  shear  deformation  as  well  as  fiber  compressibility  are  presented  in 
Sections  3  and  4,  respectively.  Section  5  talks  about  the  development  of  constitutive  matrices 
for  composite  laminates.  The  underlying  variational  framework  and  the  ensuing  finite  element 
formulation  of  the  problem  are  presented  in  Section  6.  Section  7  presents  some  numerical 
examples  to  demonstrate  the  range  of  applicability  and  accuracy  of  the  proposed  theory,  and 
conclusions  are  drawn  in  Section  8. 


2.  Assumptions  in  the  Layer- wise  Shear  Deformable  Shell  Theory 
1.  The  domain  is  of  the  following  special  form 


a=  \{x,y,z) 


zl  i 

T’  2 


(0 


where  is  the  area  of  the  reference  surface  for  layer  I,  is  the  thickness  of  layer  I  and 
T  is  the  total  thickness  of  the  composite  shell. 


2.  The  displacement  field  is  assumed  to  take  the  following  form 

U^\x,y,z)  =  (x,y,z)  +  (x,y)  a  =  1,2 


(2) 
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where  Ua\x,y,z)  are  the  in-plane  translations  in  layer  I,  9a\x,y)  are  the  out  of  plane 
director  rotations  for  the  reference  surface  associated  with  layer  I,  and  is  a  function 
that  establishes  the  position  of  a  point  with  respect  to  the  reference  surface  in  layer  1. 
The  second  term  on  the  right  hand  side  of  (2)  can  be  viewed  as  an  enhancement  to  the  in¬ 
plane  components  of  the  displacement  field  in  layer  I  with  respect  to  the  reference  surface 
associated  with  the  layer.  Consequently  this  enhanced  field  becomes  identically  zero  on 
the  reference  surfaces. 

3.  The  displacement  field  in  the  thickness  direction  is  assumed  to  be  a  function  of  z,  i.e., 

U^\x,y,z)  =  v!'^\x,y,z)  (3) 

thereby  producing  through  the  thickness  strains  which  result  in  thickness  variation  in  the 
shell.  A  consequence  of  the  above  relaxation  is  that  5^  0  i.e.,  we  do  not  invoke  the 
plane  stress  hypothesis. 

Remarks: 

1.  The  continuum  requirement  on  the  displacement  field  at  the  interface  between  and 
(Z  +  bounded  layers  necessitates  the  satisfaction  of  the  contact  conditions;  Ua^  = 

(for  q:  =  1,2),  and  These  conditions  are  inherently  satisfied  via  the 

definition  of  the  assumed  displacement  field  as  defined  in  (2)  and  (3). 

2.  The  inter-laminar  stress  continuity  requires  the  satisfaction  of  the  following  conditions  at 

the  interface  of  Z*^  and  (Z  + 1)®*  layers;  =  t^z^^  (for  a  =  1, 2),  and  where 

denotes  the  inter-laminar  stresses  for  layer  Z.  It  is  important  to  note  that  the  assmned 
displacement  field  (2)  results  in  independent  non-normal  cross-sectional  rotations  in  each 
finite  element  layer  which  is  in  accordance  with  the  MindUn  kinematic  assumption.  A  con¬ 
sequence  is  the  independent  shear  deformation  of  the  director  in  each  layer,  thus  allowing 
the  warping  of  the  composite  cross-section.  Therefore  the  resulting  strain  field  is  discon¬ 
tinuous  across  different  material  sets  and  creates  the  provision  of  stress  continuity  across 
the  interfaces.  Consequently,  the  traction  continuity  equations  are  inherently  satisfied. 
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3.  Kinematic  Description  of  Multi-layered  Shells 

Figure  1  shows  a  composite  laminate  with  ‘m’  number  of  material  layers.  (For  ease  of 
presentation,  in  Fig.  1,  m  =  3).  Consequently,  the  total  thickness  T  is  divided  into  ‘m’ 
elements  through  the  thickness.  Each  layer  of  elements  is  associated  with  a  reference  surface 
which,  for  ease  of  discussion,  is  assumed  to  be  coincident  with  the  lower  surface  of  that 
layer.  It  is  important  to  note  that  the  director  rotation  9^  and  the  slope  us,a  (^3  is  the 
displacement  field  1x3  projected  onto  the  reference  surface  for  the  layer)  are  not  necessarily 
equal  and  thus  transverse  shear  strains  are  accommodated.  This  is  to  be  contrasted  with 
the  classical  lamination  theories  in  which  6a  =  ^s.a-  Within  each  layer  the  director  rotates 
by  9^a\  generating  shear  strain  7^3.  Consequently,  in  the  deformed  configuration,  node  2 
(see  Fig.  1)  moves  to  location  2'.  Now,  the  director  in  the  second  layer  rotates  by  an  angle 
9a'^^\  generating  shear  strain  This  kinematics  is  repeated  in  successive  layers  and, 

due  to  the  continuity  of  the  displacement  field  in  the  thickness  direction,  we  obtain  the  new 
locations  of  nodal  points  (that  constitute  the  layerwise  directors)  as  l',  2',  S'  and  4'.  This  new 
location  of  points  produces  a  higher-order  variation  of  the  in-plane  displacement  field  through 
the  thickness.  Accordingly,  arbitrary  polynomial  expressions  are  not  required  to  model  the 
higher-order  in-plane  variation  of  displacement  through  the  thickness  of  the  composite. 

3.1  Kinematics  of  an  individual  layer  in  the  context  of  the  finite  element  method 
The  displacement  field  of  each  layer  I  in  the  shell  is  defined  by  the  following  relations 

n.  C)  =  (?.  n)  +  (?,  n,  c)  (4) 

^cn 

a=l 

Tien 

uPiO  (6) 

a=l 

(C)  =  ^  (0  (no  sum) 


(7) 
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Fig.  1.  Shell  kinematics.  Variation  of  transverse  shear  strains  through  the  thickness. 

4'^  =  5  (1  -  ^  +  I  (1  +  0  (S) 

where  and  are  the  translations  of  the  upper  and  lower  surfaces  of  layer  Z,  u^^'>  is  the 
displacement  of  a  point  on  the  reference  surface  of  layer  Z,  is  the  ‘director  displacement 
for  layer  Z,  Na  represents  two-dimensional  shape  function  associated  with  node  ‘a  ,  Uen 
the  number  of  element  nodes,  and  is  the  thickness  function  defined  as 

Zf  (C)  =  5  (1  +  0  -  5  (1  -  C)  Z«'  (9) 

where  ”  C)  2“’’  =  5  (1  ZW  =  xi''*  -  x<'>- 

and  =  WZ^^m  is  the  thickness  of  layer  Z.  The  parameter  C  €  [-1,+1]  defines  the  location 
of  the  reference  surface.  For  example  if  C  ~  "fjOj+l  (respectively),  then  the  reference  surface 
is  taken  to  be  the  bottom,  middle,  and  top  (respectively)  of  the  shell  layer  Z. 

The  nodal  director  displacement  vector  is  constructed  so  that  the  director  can  rotate  and 
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stretch  as  shown  by  an  additive  decomposition; 

u(i)^up+uii)  (10) 

The  through  thickness  stretch  component  can  be  expressed  as 

=  (11) 

where  and  ui‘f  are  the  translations  in  the  thickness  direction  of  the  upper  and  lower 
surfaces  belonging  to  layer  1.  The  out  of  plane  rotation  component  of  the  director  is  con- 
structed  such  that  it  may  rotate,  viz., 

t><«  =  «eir  -  (12) 

The  quantities  and  0^2  represent  the  rotations  of  the  director  about  the  basis  vectors 
and  for  shell  layer  I,  respectively. 

Remark:  We  have  used  the  right-hand-rule  sign  convention  for  rotations. 

4.  Geometric  Description  of  Multi-layered  Shells 

Figure  2  shows  a  typical  configuration  of  a  doubly  curved  composite  shell.  It  can  be  made 
of  numerous  plies  with  variable  material  properties,  ply  orientations  and  thicknesses  and  can 
be  stacked  together  in  any  arbitrary  sequence.  Figure  3  shows  a  schematic  diagram  of  the 
geometry  of  a  typical  laminate  of  the  composite  shell  (shown  by  the  shaded  region  in  Fig.  2) 
and  can  be  defined  by  the  following  relations 

X<'> C)  =  («.»))  +  (« .  0  (13) 

f^en 

a=l 

Tien 

X<‘)  «.’!.<)  =  El^«(«.'')l'”(«)  ^“) 

a=l 

(C)  =  4'^  (C)  (no  sum) 


(16) 
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where  is  the  position  vector  of  a  generic  point  in  layer  I,  is  the  position  vector  of 
a  point  on  the  reference  surface  for  layer  I,  is  the  position  vector  of  a  generic  point 
relative  to  x^^^  that  defines  the  director  through  the  point  for  layer  /,  (in  computational  shell 
literature,  X^*^  is  referred  to  as  fiber  direction),  xi  ^  is  a  unit  vector  emanating 

from  node  ‘a’  in  the  director  direction,  and  is  the  thickness  function  associated  with  node 
‘a’  as  defined  in  (9).  This  function  is  defined  by  the  location  of  the  reference  surface,  xi*^  is 
the  position  vector  of  nodal  point  ‘a’  in  layer  I  and  is  defined  as 

*1"  =  5  (1  -  <)  4'’"  +  5  (1  +  C) 


(17) 
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5.  Constitutive  Relations 

Most  composites  are  made  of  a  repeated  sequence  of  laminates  that  have  the  same  material 
properties  but  are  oriented  at  +9  and  -6  degrees  with  respect  to  the  extension  direction.  Let 
represent  the  linear  constitutive  matrix  for  an  individual  laminate  with  regard  to  its 
mutually  perpendicular  planes  of  elastic  symmetry.  This  constitutive  matrix  for  a  laminate 
can  be  projected  from  its  mutually  perpendicular  planes  of  elastic  symmetry  onto  the  global 
composite  coordinate  system  via  a  transformation  matrix  which  is  a  function  of  the  angle 

6.  The  transformed  constitutive  matrix  is  obtained  as 

(18) 


6.  The  VariationcJ  Framework 

We  start  with  the  Hu-Washizu  variational  principle  that  uses  the  displacement,  strain 
and  stress  fields  as  independent  variables  and  thus  provides  the  most  general  framework  to 
formulate  the  problem  and  derive  its  weak  form.  We  denote  the  strain  field  by  e  and  its 
variation  by  a.  Similarly,  we  denote  the  stress  field  by  cr  and  its  variation  by  /3.  Since  no 
boundary  conditions  are  specified  on  the  strain  or  stress  fields,  the  solution  space  and  the 
space  of  admissible  variations  coincide  for  both  strains  and  stresses. 

where  Usd  denotes  the  number  of  spatial  dimentions  and  for  the  present  case  =  3. 

The  Hu-Washizu  formulation  can  be  obtained  by  finding  the  stationary  points  of  the 
functional 

nHw(^)  '**)  •”  77  [  S’  C  :  s  dO,  f  cr  •  [V*tt  e]  dCl  —  next('*^)  (1^) 

2  Jn  Jo. 

where  the  first  term  is  the  stored  energy  function,  while  the  second  term  is  a  Lagrange 
multiplier  cr  enforcement  of  the  strain-displacement  relations.  The  assumed  displacement 
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field  (2)-(3)  results  in  an  enhanced  strain  field  because  layerwise  rotations  result  in  enhancing 

the  shear  strain  components.  We  write  the  strain  field  as  e  =  +  e,  where  V^-u  is  the 

symmetric  compatible  strain  and  i  is  the  enhanced  strain.  The  corresponding  strain  variation 
is  a  =  V^u;  d-o:,  where  is  the  compatible  and  d  is  the  enhanced  strain  variation.  Since 
no  boundary  conditions  are  prescribed  on  the  augmented  strains  as  well,  both  spaces  i.e.,  e 
and  d  conincide. 

S:= 

We  can  now  substitute  this  enhanced  strain  field  in  (19)  to  obtain  an  enhanced  strain 
version  of  the  Hu-Washizu  formulation. 

nHw(£(0),o-,u)  :=  i  f  {V^u  +  i{e)):C  :iV^u  +  i{0))dn-  [  a- ■  i dCl  -  Ue^t{u)  (20) 
2  jQ  Ja 

Following  Simo  and  Hughes  [22],  we  assume  a  priori,  an  orthogonality  constraint  on  the  stress 
field  or  with  respect  to  the  enhanced  strain  field  i{0).  This  causes  the  Lagrange  multiplier, 
i.e.,  (T,  to  drop  out  from  the  formulation,  thereby  leading  to  a  modified  displacement  model 

n(n,  0)  :=  ^  /  (V^U  +  i{0))  :  C  :  {V^u  +  i{0))  d^  -  next(«,  0)  (21) 

2  Ja 

The  spaces  relevant  to  the  modified  displacement  formulation  are 

S  =  {{u,0)  eH\n),  {u,0):  ^l-^'Jl^’'^;u  =  uonTu,e  =  0onTg}  (22) 

V  =  {(u;,'0)  e  (m,-0)  :  Q  (23) 

where  S  is  the  space  of  trial  displacements  and  trial  rotations,  and  V  is  the  associated  space 
of  weighting  functions,  respectively.  H^{Q)  denotes  the  space  of  square-integrable  functions 
along  with  their  generalized  derivatives  defined  over  Cl,  and  HQ{d)  is  the  subset  of  H^{Cl) 
whose  members  satisfy  zero  essential  boundary  conditions. 

For  the  case  of  small  strains,  assuming  the  composite  laminate  to  occupy  a  region 
in  minimization  of  (21)  leads  to  the  following  weak  form:  Given  f  :  Cl 
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g  :Tg^  h:Yh-*  find  {^,0}  €  S  such  that  for  all  {u,'tp}  €  V 


n(u,  6]  w,  V>)  =  /  (V®a;  +  a{9))  :  C  :  (V^u  +  i{9))  dO. 
Jq 


/  uj'fdQ  -  V  /  uj^hdT 

(24) 


where  V^(')  is  the  symmetric  gradient  operator,  f  denotes  the  body  force  vector  defined  over 
Cl,  h  represents  the  prescribed  tractions  (normal  and  shear)  over  boundary  r^,  and  g  are 
the  prescribed  displacements  and  rotations  defined  over  boundary  Fg.  Here  Fg  =  U  F^ 
with  F„  being  the  boundary  with  prescribed  displacements  and  F^  being  the  boundary  with 


prescribed  rotations. 

Now  consider  a  shell  consisting  of  ^Tn,’  generally  directional  but  perfectly  bonded  layers, 
and  reference  surface  of  each  layer  discretized  into  ''Uumei  finite  elements.  Let  ^  ^  — 
^e(o  X  [-jt/2,  fi(')  =  Ue=T'  and  n  =  element  domain 

in  layer  1.  We  can  write  the  principle  of  virtual  work  (24)  for  the  layered  but  perfectly  bonded 
anisotropic  composite  laminates  as 


(Vw  +  ct{9)f'>  :  :  (V"u  +  i{9)f'>  dQ 

(t) .  f{i)dCl  -  ^ 

1=1 


(25) 


Let  7  =  Vua  —  0  be  the  shear  strain  vector,  where  us  is  the  displacement  field  projected 
onto  the  reference  surface  for  layer  1.  Similarly,  let  k  =  be  the  curvature,  which  is 
also  defined  layer  wise.  We  can  write  the  energy  functional  (25)  for  the  layered  anisotropic 

laminates  as 


where  and  are  the  virtual  shear  strains  and  curvatures,  respectively.  7^')  =  - 

and  while  u>3  is  the  virtual  displacement  field  projected  onto  the  reference 
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surface  for  layer  1.  The  assumed  displacement  field  in  (2)  and  (3)  leads  to  a  modified  functional 
n,  defined  as 

n(n,0;a;,^)  =  f'|  /  /  7^')  :  C?)  :  dA 

+  J  dA-  J  ^  J  (27) 

where  (•)  dQ  is  a  volume  integral  and  (•)  dA  represents  an  area  integral. 

We  add  the  following  boundary  conditions  to  the  energy  functional  (27)  as 


on 

r.. 

a  =  l,2 

(28) 

U3  =  U3 

on 

Fu3 

(29) 

Oa  =  Occ 

on 

F^. 

a  =  1,2 

(30) 

where  Tu  —  Pit,,  U  represents  the  boundary  where  translation  boundary  conditions  are 
applied,  and  Fg^  corresponds  to  the  boundary  with  prescribed  rotations.  Ua,U3,sjid0a  rep¬ 
resent  the  prescribed  displacements  and  rotations,  respectively.  In  addition,  we  prescribe  the 
tractions  as 


T  n  ^  T 

on 

Fr 

(31) 

qn  =  q 

on 

r. 

(32) 

ra  n  ^  Th 

on 

Tm 

(33) 

T  and  q  represent  the  normal  and  shear  tractions,  respectively,  while  rh  represents  the  pre¬ 
scribed  moments.  This  leads  to  a  decomposition  of  F/,  as  F/,  =  F^LlFgUFm-  FrUFg  represents 
the  boundary  where  traction  boundary  conditions  are  applied,  and  F^  corresponds  to  the 
portion  of  boundary  with  applied  moments. 

The  strain  vectors  for  the  proposed  composite  element  are 


;  7  = 
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whers  Sj'y, /c  dind  s  r6pr6seiit  the  in-plciiie,  shear,  bending  and  through- thickness  strains, 

respectively.  We  can  combine  the  in-plane  strains  with  through-thickness  strains  to  yield 

.  =  {4'>4'4'24>)^  (35) 


The  stresses  corresponding  to  the  strains  are 


-(0 
^22 

-(0 
^33 

-(0 
^12 


q  = 


m 


m 


m 


m 


(0 

’ll 

(0 

'22 

(0 

12 


(36) 


where  r  represents  the  combined  in-plane  and  the  through-thickness  stresses,  and  q  and  m 
represent  the  shear  and  bending  stresses,  respectively. 

The  element  stiflEness  matrices  and  force  vectors  emanating  from  the  weak  form,  that  are 
obtained  via  Galerkin  approximation  and  introduction  of  element  shape  functions  can  be 
written  as 

m 

K  -  (37) 

i=i 

where  A"””*'  is  the  finite  element  assembly  operator.  The  element  stiffness  matrix  can  ex¬ 
plicitly  be  written  as 

=  f  dQ+  f  7^'^  :  :  7^*^  dA+  [  dA 

Ja(i)  ' - ^ - '  Jaw' - p' - '  Jaw'  -  ' 

modified  membrane  shear  bending 

(38) 

where  Cm  ,  C^^  and  C)^^  are  the  modified-membrane,  shear  and  bending  constitutive  matrices 
for  layer  1,  respectively.  The  right  hand  side  force  vector  can  be  written  as 


/* 


-h  •  hS^^dr 

i=i 


(39) 


where  the  first  term  is  to  the  body  force  vector  and  the  second  term  is  the  traction  vector. 
Remark:  We  have  used  selective  reduced  integration  concepts  for  the  evaluation  of  the  element 


stiffness  matrices. 
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7.  Numericed  Examples 

7.1  Free  edge  boundary- value  problem.  [45,-45]s 

The  first  numerical  simulation  is  a  prismatic  symmetric  laminate  having  traction-free 
edges  at  y  =  ±a  and  surfaces  Z  =  ±h,  and  is  subjected  to  strain  increment  only  on  its  ends 
at  X  =  constant.  Each  layer  is  composed  of  uni-directional  fiber-reinforced  material  such  that 
the  fiber  direction  is  defined  by  its  angle  9  with  respect  to  the  X-axis.  The  elastic  properties 
used  in  this  simulation  are  taken  from  Pagano  [12],  p.  4.  The  laminate  consists  of  four  uni¬ 
directional  fibrous  composite  layers,  two  with  axis  of  elastic  symmetry  at  -1-45  and  two  at  —45 
to  the  longitudinal  laminate"  axis.  (See  Fig.  4.)  1%  strain  in  opposite  directions  is  applied  at 
X  =  ±L,  while  it  is  restrained  to  move  in  the  axial,  lateral  and  thickness  directions  at  X  =  0. 
The  physical  dimensions  for  the  numerical  simulation  are  X  =  60,  Y  =  20,  Z  =  2.5,  with  12 
elements  in  X  direction,  40  elements  in  the  Y  direction  and  4  elements  through  the  thickness. 
For  each  finite  element  layer  through  the  thickness,  the  reference  surface  is  assumed  to  be 
coincident  with  the  bottom  surface  of  the  layer. 

Figure  5  shows  the  axial  displacement  distribution  at  the  top  free  surface  of  the  section 
cut  at  X  =0,  and  the  results  are  compared  with  Moires  et  al.  and  Pagano  et  al.  [12]. 
Figure  6  presents  the  inter-laminar  shear  stress  at  material  interface  and  their  corresponding 
values  from  Reddy  [20].  Figure  7  shows  the  major  stress  components  cru,  cr  12,  <ti3,  and  their 
corresponding  values  from  the  numerical  simulations  of  Pipes  et  al.  [16].  In  order  to  complete 
the  discussion,  we  have  also  presented  the  minor  stress  components  in  Fig.  8. 
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Figure  4.  Laminate  geometry  and  the  coordinate  system. 


Normalized  width 

Figure  5.  Axial  displacement  distribution  at  the  top  free  surface. 

7.2  Free  edge  boundary-value  problem  with  a  cylindrical  hole.  [45,-45]s 

The  composite  laminate  considered  in  this  simulation  has  the  same  physical  dimensions, 
material  properties  and  loading  conditions  as  in  the  previous  case.  However  the  present 
laminate  has  a  unit  diameter  cylindrical  hole  at  (0,0,0)  with  its  axis  coincident  with  the 
Z  axis  (see  Fig.  9).  The  ratio  of  the  diameter  of  the  hole  to  the  width  of  the  laminate  is 
0.1.  Traction-free  boundary  conditions  are  applied  on  the  cylindrical  surface  of  the  hole.  The 
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Normalized  width 


Figure  6.  Interlaminar  shear  stress  at  material  interface. 


Normalized  width 

Figure  7.  Major  stress  components  in  the  top  +45  laminate. 

laminate  is  constrained  to  move  in  the  axial,  lateral  and  transverse  directions  by  appropriately 
constraining  the  nodes  along  the  S3nnmetry  lines  at  X  =  0. 

Figure  10  presents  the  major  stress  components  <rii,  cri2,  o’la  at  a  section  cut  at  A"  =  0 
in  the  top  layer  with  +45  degrees  ply  orientation.  These  results  have  been  compared  with 
complete  3D  anisotropic  calculations  done  with  a  mesh  which  is  twice  as  refined  as  the  present 
mesh.  In  the  region  away  from  the  free  edge  boundary  and  the  traction  free  hole,  the  ratio  of 
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Normalized  width 

Figure  8.  Minor  stress  components  in  the  top  +45  laminate. 

(711  and  a  12  agrees  closely  with  that  of  the  preceding  numerical  simulation.  Figure  11  shows 
the  minor  stress  components  which  also  show  a  good  agreement  with  the  3D  anisotropic 
simulation. 


Figure  9.  Laminate  geometry  with  a  cylindrical  hole. 
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Normalized  width 


Figure  10.  Major  stress  components  in  the  top  +45  laminate. 


Normalized  width 

Figure  11.  Minor  stress  components  in  the  top  +45  laminate. 

8.  Conclusions 

In  this  paper  we  have  presented  a  finite  element  formulation  which  is  suitable  for  the 
analysis  of  multi-layered/multi-director  and  shear-deformable  composite  laminates.  The  ge¬ 
ometric  description  employed  for  the  composite  shells  finds  its  roots  in  the  so-called  degener¬ 
ated  shell  approach.  A  set  of  kinematic  hypothesis  is  introduced  to  accommodate  the  effects 
of  transverse  warping  of  the  cross  section  due  to  shear  deformation  along  with  fiber  com- 
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pressibility  that  results  from  transverse  normal  stresses.  The  kinematics  are  individually  and 

independently  represented  for  each  layer  and  are  coupled  via  the  conditions  of  displacement 
continuity  between  contiguous  layers.  The  rotation  field  is  however  layer-wise  continuous  and 
is  assumed  discontinuous  across  the  layers.  Transversal  warping  of  the  composite  cross-section 
is  described  by  individual  layer  directors  which  simultaneously  rotate  and  stretch.  This  re¬ 
sults  in  discontinuous  strain  fields  across  the  different  material  sets,  thereby  creating  the 
provision  of  stress  continuity  across  the  material  interfaces.  Since  the  formulation  resolves 
all  strains,  all  stresses  are  computed  through  the  three-dimensional  constitutive  equations 
and  the  usual  ‘zero  normal  stress’  shell  hypothesis  is  not  employed.  The  variational  equation 
contains  only  the  first  derivatives  of  displacement  and  rotation  fields  that  require  just  the  C® 
continuity  of  finite  element  functions  [7|.  In  displacement  formulation  of  plates  and  shells, 
special  attention  needs  to  be  given  to  transverse  shear  and  membrane  terms  to  prevent  the 
occurrence  of  Tnesh  locking.  We  have  employed  the  selective/reduced  integration  technique  to 
avoid  this  problem.  Numerical  results  are  presented  that  demonstrate  the  good  performance 
of  the  proposed  formulation. 
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